Tidy data

Set-up

knitr::opts_chunk$set(
    echo = TRUE,
    error = FALSE,
    message = FALSE,
    warning = FALSE
)
library(afex)
library(tidyverse)
library(kableExtra)
library(knitr)
library(plotly)
library(broom)
library(wesanderson)
library(robust)
library(ggparl)
library(mediation)

set.seed(2022)

Read data

# reads excel data sheet

data <- read.csv("../data/study4.csv") # reads data
data_sat <- read.csv("../data/raw/stage3_satisfaction.csv") # reads in satisfaction data



# summary(data) # explore the data types, missing data and typos

data <- data %>%
  mutate(gender = ifelse(gender == "m", "m", "f")) # corrects an issuewith extra space in "f "

data$gender <- as.factor(data$gender)
data$pp_number <- as.factor(data$pp_number)
data$condition <-as.factor(data$condition)


# mood data
mood1 <- read.csv("../data/raw/stage1_mood.csv") %>%
  dplyr::select(subject, stimulusitem2, response)%>%
  rename(response1=response, feeling = stimulusitem2)%>%
  mutate(stage = 1)%>%
  filter(subject != 6666)

mood2 <- read.csv("../data/raw/stage2_mood.csv") %>%
  dplyr::select(subject, stimulusitem2,  response)%>%
  rename(response2=response,  feeling = stimulusitem2)%>%
  mutate(stage = 2)%>%
  filter(subject != 6666)

mood3 <- read.csv("../data/raw/stage3_mood.csv") %>%
  dplyr::select(subject, stimulusitem2, response) %>%
  rename(response3=response,  feeling = stimulusitem2)%>%
  mutate(stage = 3)%>%
  filter(subject != 6666)

mood_dat <- inner_join( mood1, mood2,
                        by = c("subject", "feeling"), 
                        keep = FALSE)%>%
            inner_join(., mood3,
                       by = c("subject", "feeling"),
                       keep = FALSE)%>%
            dplyr::select(- starts_with("stage")) %>%
            pivot_longer(cols = response1:response3,
                         names_to ="stage",
                         values_to = "rating")%>%
            pivot_wider(names_from = "feeling",
                        values_from ="rating" )%>%
            mutate(alert = (alert +strong + `clear-headed` + `well-coordinated` + energetic + `quick-witted` + attentive + proficient + interested )/9,
                   content =(contented + tranquil + amicable + gregarious)/4,
                   calm =(calm + relaxed)/2)




#`clear-headed` + `well-coordinated` + energetic + `quick-witted` + attentive + proficient + interested

# qualtrics data 
      #mainly info about drinking

data_q<- read.csv("../data/other/study_4Q.csv")%>%
  dplyr::select(Q21, gender, age, Q15, Q16)%>%
  rename(subject = Q21, 
         units = Q15, 
         units_beer = Q16)%>%
  mutate(subject = as.numeric(subject),
         condition = subject %/% 1000)%>%
  filter(!is.na(subject))%>%
  mutate(subject = as.factor(subject), 
         gender=as.factor(gender), 
         age = as.numeric(age), 
         units = as.numeric(units), 
         units_beer = as.numeric(units_beer))%>%
  filter( subject != "6666")


# expectations & perception of taste/flavour

taste_dat <- read.csv("../data/raw/stage1_perception.csv")%>%
  dplyr::select(subject, trialcode, response)%>%
  rename(Epercept = trialcode, expected = response)

taste_dat$Epercept<- as.factor(taste_dat$Epercept)
  

expect_dat<- read.csv("../data/raw/stage1_expectations.csv")%>%
  dplyr::select(subject, trialcode, response, stimulusitem3)%>%
  filter(stimulusitem3 == "extremely")%>%
  dplyr::select(- stimulusitem3)%>%
  rename(percept = trialcode, perceived = response)

expect_dat$percept<- as.factor(expect_dat$percept)


# satisfaction data
  #merge data_sat and expect_dat

satis_datP <- data_sat %>% 
  dplyr::select(subject, trialcode, response)%>%
  dplyr::filter(trialcode == "VAS")

satis_datE <- expect_dat %>%
  dplyr::filter(percept == "VAS_satisE_VAS")


wtp_dat <- data_sat %>%
  filter(trialcode == "VAS_wtp")%>%
  dplyr::select(subject, response)


# ITT data

itt_s1 <-read.csv("../data/summaries/stage1_ITT.csv")%>%
  dplyr::select(script.subjectid, expressions.propcorrect_stim1:expressions.propcorrect_stim13)%>%
  mutate(stage1 = 1)%>%
  rename(subject = script.subjectid)%>%
  rename_at(vars(matches("^expressions")), ~ str_remove(.,"^expressions.propcorrect_")) %>%
  pivot_longer(cols = starts_with("stim") , 
               names_to = "stimulus" ,
               values_to = "correct")


itt_s2 <-read.csv("../data/summaries/stage2_ITT.csv")%>%
  dplyr::select(script.subjectid, expressions.propcorrect_stim1:expressions.propcorrect_stim13)%>%
  mutate(stage2 = 2)%>%
  rename(subject = script.subjectid)%>%
  rename_at(vars(matches("^expressions")), ~ str_remove(.,"^expressions.propcorrect_"))%>%
  pivot_longer(cols = starts_with("stim") , 
               names_to = "stimulus" ,
               values_to = "correct")

itt_s3 <-read.csv("../data/summaries/stage3_ITT.csv")%>%
  dplyr::select(script.subjectid, expressions.propcorrect_stim1:expressions.propcorrect_stim13)%>%
  mutate(stage3 = 3)%>%
  rename(subject = script.subjectid)%>%
  rename_at(vars(matches("^expressions")), ~ str_remove(.,"^expressions.propcorrect_"))%>%
  pivot_longer(cols = starts_with("stim") , 
               names_to = "stimulus" ,
               values_to = "correct")

itt_dat <- inner_join(itt_s1, itt_s2,
                      by = c("subject","stimulus"),
                      keep = FALSE,
                      suffix = c(".1", ".2"))%>%
          inner_join(., itt_s3,
                     by = c("subject", "stimulus"),
                     keep = FALSE) %>%
  pivot_longer(cols = starts_with("correct") ,
               names_to = "stage",
               values_to = "correct")%>%
  dplyr::select(-c( stage1, stage2, stage3))%>%
  mutate( stage = case_when(stage == "correct.1" ~ 'stage 1',
                            stage == "correct.2" ~ "stage 2",
                            stage == "correct" ~ "stage 3"))%>%
   mutate(condition = subject %/% 1000) %>%
  mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
  mutate(focus = case_when(condition == 1 ~ 'control',
                           condition == 2 ~ 'control',
                           condition == 3 ~ 'hedonic',
                           condition == 4 ~ 'hedonic',
                           condition == 5 ~ 'utilitarian',
                           condition == 6 ~ 'utilitarian'))

Data joins

DA <- data %>%
  dplyr::select(pp_number, BMI)%>%
  rename(subject = pp_number)

DB <- data_sat %>%
  dplyr::select(subject, response, trialcode ) %>%
  mutate(row = row_number(), subject = as.factor(subject)) %>%
  pivot_wider(names_from = trialcode, values_from = response) %>%
  rename( wtp = VAS_wtp)%>%
  dplyr::select(- row)


D1 <- inner_join( DA,DB,
                  by = "subject",
                  keep = FALSE
                ) %>%
    pivot_longer(cols =c(VAS, wtp) , 
                 names_to = "measure", 
                 values_to = "rating") %>%
    inner_join(., data_q,
             by = "subject", 
             keep = FALSE)%>%
  filter(!is.na(rating))
                                            # long dataset incl,. info about satisfaction and willingness to pay by condition, and alcohol/ non-alcohol
# need to add focus condition


#Joining mood data

D2 <- inner_join(mood1, mood2, 
                 by = c("subject", "feeling"),
                 keep = FALSE) %>%
  inner_join(., mood3,
             by = c("subject","feeling"),
             keep = FALSE)%>%
            dplyr::select(- starts_with("stage")) %>%
  mutate(condition = subject  %/% 1000) %>%
  pivot_longer(cols = starts_with("response"),
   names_to = "stage",
   names_prefix = "response",
   values_to = "rating")
  

D3_dat <- inner_join(expect_dat, taste_dat,
                     by = "subject",
                     keep = FALSE)%>%
  mutate(condition = subject %/% 1000) %>%
  pivot_wider(names_from = percept,
    values_from = perceived) %>%
  pivot_wider(names_from = Epercept,
              values_from = expected)%>%
  rename_at(vars(matches("^VAS\\_")), ~ str_remove(.,"^VAS\\_"))%>%
  rename_at(vars(matches("_VAS")), ~ str_remove(.,"_VAS")) %>%
  mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
  mutate(focus = case_when(condition == 1 ~ 'control',
                           condition == 2 ~ 'control',
                           condition == 3 ~ 'hedonic',
                           condition == 4 ~ 'hedonic',
                           condition == 5 ~ 'utilitarian',
                           condition == 6 ~ 'utilitarian'))
  


D3_dat$alcohol<- as.factor(D3_dat$alcohol)
D3_dat$focus<- as.factor(D3_dat$focus)
D3_dat$condition <- as.factor(D3_dat$condition)


satis_dat <- inner_join(satis_datE, satis_datP,
                        by = "subject",
                        keep = FALSE)%>%
              dplyr::select(subject, perceived, response,)%>%
            rename(satisE = perceived, satis = response)%>%
              inner_join(., wtp_dat,
                         by = "subject",
                         keep = FALSE)%>%
            rename(wtp = response) %>%
  mutate(condition = subject %/% 1000) %>%
  mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
  mutate(focus = case_when(condition == 1 ~ 'control',
                           condition == 2 ~ 'control',
                           condition == 3 ~ 'hedonic',
                           condition == 4 ~ 'hedonic',
                           condition == 5 ~ 'utilitarian',
                           condition == 6 ~ 'utilitarian'))

Descriptive analyses

Differences between conditions

# gender
total_n <- nrow(data)
males <- data_q%>%filter(gender == "male") %>% nrow()
females <- data_q%>%filter(gender == "female") %>% nrow()
nonbinary <- data_q %>% filter(gender== "non-binary" )%>% nrow()



# by condition
gender_tib <- D1 %>% 
  filter(measure == "wtp") %>% # avoids duplicating data
  group_by( condition) %>%
  summarise(n=n(),
            BMI = mean(BMI, na.rm = TRUE), 
            age = mean(age, na.rm = TRUE)) %>% 
  kable()

Gender

xtabs(~ gender + condition , data = data_q,drop.unused.levels = TRUE)%>%
  kable(digits = 2, caption = "A contingency table with observed gender values/counts across conditions") %>%
  kable_styling(full_width = TRUE, 
                position = "left")
A contingency table with observed gender values/counts across conditions
1 2 3 4 5 6
female 24 21 22 25 30 23
male 6 8 7 6 6 8
non-binary 2 1 1 1 0 0
chi_test <- chisq.test(data_q$gender, data_q$condition, simulate.p.value = TRUE)

chi_test$expected %>%
  kable(digits = 1, caption = "A contingency table of expected gender values/counts across conditions")%>%
  kable_styling(full_width = TRUE,
                position = "left")
A contingency table of expected gender values/counts across conditions
1 2 3 4 5 6
female 24.3 22.8 22.8 24.3 27.3 23.5
male 6.9 6.4 6.4 6.9 7.7 6.7
non-binary 0.8 0.8 0.8 0.8 0.9 0.8
#  χ2 test requires all expected frequencies to be greater than five --> simulate.p.value avoids the problem

Altogether 189 participants took part in the study, as is common for psychology studies, most of these were women (145 self-identified as female, 41 as male and 5 as non-binary). The gender split across the conditions was not problematic \(\chi^2\) = 5.28, p = 0.903 (simulated p values based on 2000 replicates, as expected values << 5).

# plot
## bar plot is best here
gender_pal <- c("#d62828", "#003049", "#9d69a3")
gender_cond_plot <-
  ggplot(data = data_q)+
  geom_bar(aes(x= condition, fill = gender))+
  scale_fill_manual(values = gender_pal) +
  scale_y_continuous(name="", limits=c(0, 40), breaks = seq(0,40, 5)) +
  scale_x_discrete(name = "conditions", breaks = seq(1,6,1), limits = c("1", "2", "3", "4", "5", "6"))+
  theme(panel.grid = element_blank())+
  theme_minimal()

 ggplotly(gender_cond_plot)  

Number of participants in each conditions

Age

data_desc <- D1 %>%
  filter(measure == "wtp")

age_lm <- lm(age ~ condition, data=data_desc)

age_sum <- broom::tidy(age_lm)
age_glance <- broom::glance(age_lm)

Participants’ age ranged between 18 and 57. Most participants were very young and the mean age of participants was 20.3 (sd = 4.7). The age did not significantly differ across the six conditions F(181, 1) = 0.26, p = 0.608.

age_plot <- ggplot(data= data_desc, aes(x=condition, y = age, group = condition) )+
  geom_boxplot(fill= "#6c757d", alpha = 0.5, outlier.shape = NA) +
  geom_jitter(shape=16, position=position_jitter(0.2), aes(color = gender))+
  scale_y_continuous(name="age", limits=c(18, 60), breaks = seq(0,60, 10)) +
  scale_x_discrete(name = "conditions", limits = c("1", "2", "3", "4", "5", "6"))+
  scale_color_manual(values = gender_pal) +
  theme_minimal()

age_plot2 <- ggplot(data= data_desc, aes(x=condition, y = age, group = condition) )+
  geom_boxplot(fill= "#6c757d", alpha = 0.5, outlier.shape = NA) +
  geom_jitter(shape=16, position=position_jitter(0.2), aes(color = gender))+
  scale_y_continuous(name="age", limits=c(15, 30), breaks = seq(15,30,5)) +
  scale_x_discrete(name = "conditions", limits = c("1", "2", "3", "4", "5", "6"))+
  scale_color_manual(values = gender_pal) +
  theme_minimal()

ggplotly(age_plot)

Participants’ age by condition and gender. Note range of the y-axis

ggplotly(age_plot2)

Participants’ age by condition and gender. Note range of the y-axis

BMI

bmi_lm <- lm(BMI ~ condition, data=data_desc)

bmi_glance <- broom::glance(bmi_lm)

The participants’ BMI ranged between 16.7 - 36.1. The average BMI was 23.1 (sd =3.7) and again, this did not differ across conditions F(186, 1) = 1.62, p = 0.205.

BMI_plot <- ggplot(aes(x=condition, y = BMI), data= data_desc)+
  geom_boxplot(fill= "#6c757d", alpha = 0.5) +
  geom_jitter(shape=16, position=position_jitter(0.2), aes(color = gender))+
  scale_y_continuous(name="BMI", limits=c(18, 35), breaks = seq(0,35, 10)) +
  scale_x_discrete(name = "conditions", limits = c("1", "2", "3", "4", "5", "6"))+
  scale_color_manual(values = gender_pal) +
  theme_minimal()

ggplotly(BMI_plot)

Participants’ BMI by condition and gender

Baseline mood

Visual inspection of participants’ baseline mood does not give raise to any concerns. Do I need to run 16 anovas? Does not seem necessary…)

mood_base <- D2 %>%
  filter(stage == 1)

palette <- c("#033f63","#28666e","#7c9885","#b5b682","#fedc97","#e5e1ee","#c17767","#210203","#eb9486","#82204a", "#0a2342","#2ca58d","#84bc9c","#fffdf7","#f46197","#fec601")

ggplot(data = mood_base, aes(x=condition, y = rating, group = condition, fill = feeling))+
  geom_boxplot()+
  facet_wrap("feeling")+
  scale_fill_manual(values = palette) +
  scale_y_continuous(name="", limits=c(0, 100), breaks = seq(0,100, 20)) +
  theme(panel.grid = element_blank())+
  theme_minimal()
Participants' baseline mood ratings

Participants’ baseline mood ratings

Drinking

alcohol consumption

While the average alcohol consumption was realtively low 10.6 (sd = 8.9), there was a considerable differences among participants. The lowest alcohol consumption was only 1, equivalent to a half a pint of light beer (3.5%) once a week, the highest alcohol consumption was 70 units of alcohol a week, which is equivalent to over 20 pints of strong (~5% abv) beer every week. Interestingly, 45 (23.8%) participants consumed more than the recommended 14 units a week. We also asked participants to estimate how many units of beer they consumed. This was on average 4.8 (sd = 4.7), this however varied between 0 and 30. Thefigures below show participants’ weekly alcohol consumption, weekly beer consumption and the proportion of units of alcohol that came from beer. Neither the overall alcohol consumption nor beer consumption differed significantly across the conditions, F(2, 189) = 3.19, p = 0.069 and F(2, 189) = 0.06, p = 0.804, respectively.

# plot units alcohol consumption by condition

units_plot<-ggplot(aes(x=condition, 
           y = units, 
           group = condition), 
           data= data_q)+
  geom_hline(aes(yintercept = 14, linetype = "14 units/week"),
             color = "red")+
  scale_linetype_manual(name = "max. recommended\nalcohol consumption", values = c(2, 2),
                      guide = guide_legend(override.aes = list(color = c("red")), title.theme = element_text(size = 8, face = "bold" ) ))+
  geom_boxplot(fill= "#6c757d", 
               alpha = 0.5, 
               outlier.shape = NA) +
  geom_jitter(shape=16, 
              position=position_jitter(0.2), 
              aes(color = gender))+
  scale_y_continuous(name="alcohol consumption\n(units)", 
                     limits=c(0, 75), 
                     breaks =   seq(0,75, 5)) +
  scale_x_discrete(name = "conditions", 
                   limits = c("1", "2", "3", "4", "5", "6"))+
  scale_color_manual(values = gender_pal) +
  theme_minimal()

ggplotly(units_plot)

participants’ alcohol consumption by condition

# plot units beer consumption by condition

beer_plot<-ggplot(aes(x=condition, 
           y = units_beer, 
           group = condition), 
           data= data_q)+
  geom_boxplot(fill= "#6c757d", 
               alpha = 0.5, 
               outlier.shape = NA) +
  geom_jitter(shape=16, 
              position=position_jitter(0.2), 
              aes(color = gender))+
  scale_y_continuous(name="beer consumption\n(units)", 
                     limits=c(0, 35), 
                     breaks =   seq(0,35, 5)) +
  scale_x_discrete(name = "conditions", 
                   limits = c("1", "2", "3", "4", "5", "6"))+
  scale_color_manual(values = gender_pal) +
  theme_minimal()

ggplotly(beer_plot)

participants’ beer consumption by condition

summary_beer <- data_q %>%
   pivot_longer(cols = starts_with("units"),  names_to = "al_id", values_to = "units") %>%
   mutate(al_id = ifelse(al_id =="units", "all", "beer"))%>%
   group_by(condition, al_id)%>%
   summarize(units = mean(units))%>%
  pivot_wider(names_from = al_id, values_from = units)%>%
  mutate(perc= (beer/all)*100)%>%
  pivot_longer(cols = c(all, beer), names_to = "alcohol", values_to = "units")

al_pal <- c( "#1d3557","#f4a261")
  
beer_prop <- ggplot(data=summary_beer, 
       aes(x=condition,
           y=units, 
           fill = alcohol)) +
  geom_bar(stat="identity", 
           position=position_dodge())+
  geom_text(aes(label =paste0(round(perc,1),"%") ), 
             data = summary_beer %>%filter(alcohol == "beer"), 
             show.legend = FALSE,
             size = 2, 
             nudge_x = 0.3, 
             nudge_y = 1,
             color = "#1d3557" )+
  scale_fill_manual(values = al_pal)+
  scale_y_continuous(name="alcohol consumption\n(units)", 
                     limits=c(0, 15), 
                     breaks =   seq(0,15, 5)) +
  scale_x_discrete(name = "conditions", 
                   limits = c("1", "2", "3", "4", "5", "6"))+
  labs( fill = "alcohol")+
  theme_minimal()

ggplotly(beer_prop)

% of alcohol consumption that is beer by condition

Expectations & Gustatory Perception

Bitterness

  • no direct or indirect effect of focus or alcoholcontent on perceived orexpected bitterness
bitter_dat <- D3_dat%>%
  dplyr::select(subject, alcohol, focus, condition, bitterE, bitter)
   # prep dataset

# model m A --> M
mbm <- lm(bitterE ~ focus*alcohol, data = bitter_dat)

# model A--> B
mby <-lm(bitter ~ bitterE + focus*alcohol,  data = bitter_dat)

# robustSE = TRUE to get heteroscedasticity consistent SE
# nonparametric bootstrap rather than the quasi-Bayesian Monte
# Carlo simulation for variance estimation via the boot = TRUE argument
# also entering alcohol/focus as a covariate

broom::tidy(mbm)
## # A tibble: 6 x 5
##   term                                  estimate std.error statistic  p.value
##   <chr>                                    <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)                              55.5       2.80    19.8   1.61e-47
## 2 focushedonic                              5.26      4.06     1.29  1.97e- 1
## 3 focusutilitarian                         -3.08      3.85    -0.801 4.24e- 1
## 4 alcoholnon-alcoholic                     -3.70      4.03    -0.919 3.59e- 1
## 5 focushedonic:alcoholnon-alcoholic         1.22      5.72     0.214 8.31e- 1
## 6 focusutilitarian:alcoholnon-alcoholic     6.99      5.59     1.25  2.13e- 1
med_bitter_focusH <- mediate(mbm, mby, 
                            sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "bitterE",
                            covariates = "alcohol",
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "hedonic",)

med_bitter_focusU <- mediate(mbm, mby, 
                            sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "bitterE",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "utilitarian")


summary(med_bitter_focusH)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME             1.5961       0.0275         3.82   0.044 *
## ADE              1.6530      -5.8790         9.24   0.678  
## Total Effect     3.2490      -4.2910        10.74   0.422  
## Prop. Mediated   0.4912      -3.3536         3.22   0.446  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_bitter_focusH)

summary(med_bitter_focusU)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME             0.0925      -1.5001         1.73    0.88
## ADE              3.7717      -3.0102        10.89    0.32
## Total Effect     3.8642      -3.2772        10.99    0.29
## Prop. Mediated   0.0239      -1.5064         1.58    0.86
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_bitter_focusU)

med_bitter_alcohol <- mediate(mbm, mby, 
                              sims = 1000, 
                              boot.ci.type = "perc", 
                              treat = "alcohol", 
                              mediator = "bitterE",
                              robustSE = TRUE, 
                              boot = TRUE,
                              conf.level = 0.95, 
                              covariates = "focus", 
                              control.value = "alcoholic", 
                              treat.value = "non-alcoholic")

summary(med_bitter_alcohol)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME             -0.229       -1.629         1.14    0.70
## ADE              -0.167       -6.500         6.16    0.99
## Total Effect     -0.397       -6.912         5.91    0.92
## Prop. Mediated    0.578       -1.912         2.72    0.87
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_bitter_alcohol)

bitter_tib <-bitter_dat %>%
  pivot_longer(cols = starts_with("bitter"),
               names_to = "bitter",
               values_to = "rating")%>%
  group_by(alcohol, focus,bitter)%>%
  summarise(n = n(),
            meanB = mean(rating, na.rm = TRUE),
            sdB = sd(rating, na.rm = TRUE),
            ci_lo = meanB - 1.96 * sdB/sqrt(n),
            ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
  filter(!is.na(focus))

rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("bitter", "bitterE")

bitter_plot<-ggplot(aes(x= focus, y = meanB, fill = alcohol), data = bitter_tib ) +
  facet_grid(~ bitter, 
             labeller = labeller(bitter = rating_labs))+
   geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
                 position=position_dodge(.9))+
  scale_y_continuous(name="bitterness", 
                     limits=c(0, 100), 
                     breaks =   seq(0,100, 25)) +
  scale_x_discrete(name = "focus" )+
  scale_fill_manual(values = gender_pal) +
  theme_minimal()

ggplotly(bitter_plot)

Expected and perceived ratings of bitterness errorbars indicate 95% CI (1.96 SE)

Refreshment

  • no direct or indirect effect of focus/ alcohol content on perceived refreshment
  • focus on effects of drink –> increased expectations of refreshment, but did not affect perception
refresh_dat <- D3_dat%>%
  dplyr::select(subject, alcohol, focus, condition, refreshmentE, refreshment)
   # prep dataset

# model m A --> M
mrm <- lm(refreshmentE ~ focus*alcohol, data = refresh_dat)

# model A--> B
mry <-lm(refreshment ~ refreshmentE + focus*alcohol,  data = refresh_dat)

broom::tidy(mrm)%>%
  kable(digits = 3, caption = "effect of alcohol and focus on expected refreshment")%>%
  kable_styling(full_width = TRUE)
effect of alcohol and focus on expected refreshment
term estimate std.error statistic p.value
(Intercept) 63.250 3.052 20.725 0.000
focushedonic 0.474 4.426 0.107 0.915
focusutilitarian 6.472 4.194 1.543 0.125
alcoholnon-alcoholic -2.417 4.387 -0.551 0.582
focushedonic:alcoholnon-alcoholic 2.818 6.232 0.452 0.652
focusutilitarian:alcoholnon-alcoholic -0.435 6.095 -0.071 0.943
med_refreshment_focusH <- mediate(mrm, mry, 
                                  sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "refreshmentE",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "hedonic")

summary(med_refreshment_focusH)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME              0.476       -1.180         1.99   0.518  
## ADE              -5.184      -10.226         0.31   0.064 .
## Total Effect     -4.708       -9.785         0.80   0.090 .
## Prop. Mediated   -0.101       -1.410         0.73   0.584  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_refreshment_focusH)

med_refreshment_focusU <- mediate(mrm, mry, 
                                  sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "refreshmentE",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "utilitarian")

summary(med_refreshment_focusU)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME             1.6072       0.0126         3.59   0.046 *
## ADE             -3.8968      -8.8457         1.17   0.120  
## Total Effect    -2.2895      -7.5390         3.10   0.418  
## Prop. Mediated  -0.7020     -11.5112         7.92   0.456  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_refreshment_focusU)

med_refreshment_alcohol <- mediate(mrm, mry,
                                   sims = 1000, 
                                  boot.ci.type = "perc", 
                                  treat = "alcohol", 
                                  mediator = "refreshmentE",
                                  robustSE = TRUE, 
                                  boot = TRUE,
                                  conf.level = 0.95, 
                                  covariates = "focus", 
                                  control.value = "alcoholic", 
                                  treat.value = "non-alcoholic")

summary(med_refreshment_alcohol)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME             -0.428       -1.945         0.80    0.50
## ADE              -1.600       -5.955         2.61    0.43
## Total Effect     -2.028       -6.515         2.32    0.33
## Prop. Mediated    0.211       -2.628         2.81    0.56
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_refreshment_alcohol)

refreshment_tib <-refresh_dat %>%
  pivot_longer(cols = starts_with("refreshment"),
               names_to = "refreshment",
               values_to = "rating")%>%
  group_by(alcohol, focus,refreshment)%>%
  summarise(n = n(),
            meanB = mean(rating, na.rm = TRUE),
            sdB = sd(rating, na.rm = TRUE),
            ci_lo = meanB - 1.96 * sdB/sqrt(n),
            ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
  filter(!is.na(focus))

rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("refreshment", "refreshmentE")

refreshment_plot <- ggplot(aes(x= focus, y = meanB, fill = alcohol), data = refreshment_tib ) +
  facet_grid(~ refreshment, 
             labeller = labeller(refreshment = rating_labs))+
   geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
                 position=position_dodge(.9))+
  scale_y_continuous(name="refreshment", 
                     limits=c(0, 100), 
                     breaks =   seq(0,100, 25)) +
  scale_x_discrete(name = "focus" )+
  scale_fill_manual(values = gender_pal) +
  theme_minimal()

ggplotly(refreshment_plot)

Expected and perceived ratings of Refreshment errorbars indicate 95% CI (1.96 SE)

Liking

  • no effect of focus or alcohol content on expected or perceived liking
liking_dat <- D3_dat%>%
  dplyr::select(subject, alcohol, focus, condition, likingE, liking)
   # prep dataset

# model m A --> M
mlm <- lm(likingE ~ focus*alcohol, data = liking_dat)

# model A--> B
mly <-lm(liking ~ likingE + focus*alcohol,  data = liking_dat)


broom::tidy(mlm)%>%
  kable(digits = 3, caption = "effect of alcohol and focus on expected liking")%>%
  kable_styling(full_width = TRUE)
effect of alcohol and focus on expected liking
term estimate std.error statistic p.value
(Intercept) 58.594 3.404 17.212 0.000
focushedonic -6.249 4.937 -1.266 0.207
focusutilitarian -1.455 4.679 -0.311 0.756
alcoholnon-alcoholic -5.394 4.894 -1.102 0.272
focushedonic:alcoholnon-alcoholic 2.955 6.952 0.425 0.671
focusutilitarian:alcoholnon-alcoholic -0.261 6.798 -0.038 0.969
med_liking_focusH <- mediate(mlm, mly, 
                            sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "likingE",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "hedonic")

summary(med_liking_focusH)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME             -0.745       -2.151         0.39   0.202  
## ADE              -5.091      -11.422         1.42   0.112  
## Total Effect     -5.835      -12.054         0.77   0.078 .
## Prop. Mediated    0.128       -0.465         0.89   0.256  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_liking_focusH)

med_liking_focusU <- mediate(mlm, mly, 
                            sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "likingE",
                            covariates = "alcohol",
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "utilitarian")

summary(med_liking_focusU)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME            -0.2454      -1.3588         0.92    0.67
## ADE             -3.0973      -9.3443         2.94    0.30
## Total Effect    -3.3427      -9.8839         2.99    0.27
## Prop. Mediated   0.0734      -1.3463         1.02    0.65
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_liking_focusU)

med_liking_alcohol <- mediate(mlm, mly, 
                              sims = 1000, 
                              boot.ci.type = "perc", 
                              treat = "alcohol", 
                              mediator = "likingE",
                              robustSE = TRUE, 
                              boot = TRUE,
                              conf.level = 0.95, 
                              covariates = "focus", 
                              control.value = "alcoholic", 
                              treat.value = "non-alcoholic")

summary(med_liking_alcohol)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME             -0.703       -2.170         0.11    0.12
## ADE              -3.393       -8.464         2.01    0.20
## Total Effect     -4.096       -9.395         1.24    0.11
## Prop. Mediated    0.172       -0.529         1.51    0.20
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_liking_alcohol)

liking_tib <-liking_dat %>%
  pivot_longer(cols = starts_with("liking"),
               names_to = "liking",
               values_to = "rating")%>%
  group_by(alcohol, focus,liking)%>%
  summarise(n = n(),
            meanB = mean(rating, na.rm = TRUE),
            sdB = sd(rating, na.rm = TRUE),
            ci_lo = meanB - 1.96 * sdB/sqrt(n),
            ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
  filter(!is.na(focus))

rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("liking", "likingE")

liking_plot <- ggplot(aes(x= focus, y = meanB, fill = alcohol), data = liking_tib ) +
  facet_grid(~ liking, 
             labeller = labeller(liking = rating_labs))+
   geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
                 position=position_dodge(.9))+
  scale_y_continuous(name="liking", 
                     limits=c(0, 100), 
                     breaks =   seq(0,100, 25)) +
  scale_x_discrete(name = "focus" )+
  scale_fill_manual(values = gender_pal) +
  theme_minimal()

ggplotly(liking_plot)

Expected and perceived ratings of liking errorbars indicate 95% CI (1.96 SE)

Body

  • alcoholic beer rated as having fuller body vs. non-alcoholic beer, however not mediated by expectations of body -no effect of focus on expected or perceived body
body_dat <- D3_dat%>%
  dplyr::select(subject, alcohol, focus, condition, bodyE, body)
   # prep dataset

# model m A --> M
mbom <- lm(bodyE ~ focus*alcohol, data = body_dat)

# model A--> B
mboy <-lm(body ~ bodyE + focus*alcohol,  data = body_dat)

broom::tidy(mbom)%>%
  kable(digits = 3, caption = "effect of alcohol and focus on expected body") %>%
  kable_styling(full_width = TRUE)
effect of alcohol and focus on expected body
term estimate std.error statistic p.value
(Intercept) 48.594 3.100 15.676 0.000
focushedonic 1.682 4.496 0.374 0.709
focusutilitarian 3.240 4.260 0.760 0.448
alcoholnon-alcoholic -7.727 4.456 -1.734 0.085
focushedonic:alcoholnon-alcoholic 6.857 6.330 1.083 0.280
focusutilitarian:alcoholnon-alcoholic 3.636 6.190 0.587 0.558
med_body_focusH <- mediate(mbom, mboy,
                          sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "bodyE",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "hedonic")

summary(med_body_focusH)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME              0.996       -0.152         3.36    0.14
## ADE               1.278       -5.479         7.96    0.67
## Total Effect      2.273       -4.116         9.01    0.48
## Prop. Mediated    0.438       -2.955         2.96    0.53
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_body_focusH)

med_body_focusU <- mediate(mbom, mboy,
                          sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "bodyE",
                            covariates = "alcohol",
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "utilitarian")

summary(med_body_focusU)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME             0.9918      -0.0983         3.05    0.10
## ADE             -1.6508      -8.3809         5.66    0.64
## Total Effect    -0.6590      -7.5640         6.56    0.86
## Prop. Mediated  -1.5049      -4.2345         4.24    0.91
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_body_focusU)

med_body_alcohol <- mediate(mbom, mboy,
                             sims = 1000, 
                              boot.ci.type = "perc", 
                              treat = "alcohol", 
                              mediator = "bodyE",
                              robustSE = TRUE, 
                              boot = TRUE,
                              conf.level = 0.95, 
                              covariates = "focus", 
                              control.value = "alcoholic", 
                              treat.value = "non-alcoholic")

summary(med_body_alcohol)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME            -0.8385      -2.5845         0.19   0.148  
## ADE             -6.1981     -12.1954        -0.50   0.042 *
## Total Effect    -7.0366     -13.0981        -1.26   0.016 *
## Prop. Mediated   0.1192      -0.0543         0.58   0.156  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 190 
## 
## 
## Simulations: 1000
plot(med_body_alcohol)

body_tib <-body_dat %>%
  pivot_longer(cols = starts_with("body"),
               names_to = "body",
               values_to = "rating")%>%
  group_by(alcohol, focus,body)%>%
  summarise(n = n(),
            meanB = mean(rating, na.rm = TRUE),
            sdB = sd(rating, na.rm = TRUE),
            ci_lo = meanB - 1.96 * sdB/sqrt(n),
            ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
  filter(!is.na(focus))

rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("body", "bodyE")

body_plot <-ggplot(aes(x= focus, y = meanB, fill = alcohol), data = body_tib ) +
  facet_grid(~ body, 
             labeller = labeller(body = rating_labs))+
   geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
                 position=position_dodge(.9))+
  scale_y_continuous(name="body", 
                     limits=c(0, 100), 
                     breaks =   seq(0,100, 25)) +
  scale_x_discrete(name = "focus" )+
  scale_fill_manual(values = gender_pal) +
  theme_minimal()

ggplotly(body_plot)

Expected and perceived ratings of body errorbars indicate 95% CI (1.96 SE)

Post - ingestive experience

Satisfaction

  • no effect of focus and alcohol content on expected satisfaction
  • participants rated perceived satisfaction lower when they drank non-alcoholic beer
  • effect of alcohol on perceived satisfaction not mediated by expected satisfaction
satis_tib <-satis_dat %>%
  pivot_longer(cols = starts_with("satis"),
               names_to = "satis",
               values_to = "rating")%>%
  group_by(alcohol, focus,satis)%>%
  summarise(n = n(),
            meanB = mean(rating, na.rm = TRUE),
            sdB = sd(rating, na.rm = TRUE),
            ci_lo = meanB - 1.96 * sdB/sqrt(n),
            ci_hi = meanB + 1.96 * sdB/sqrt(n))%>%
  filter(!is.na(focus))

rating_labs <- c("perceived", "expected")
names(rating_labs) <- c("satis", "satisE")

satis_plot<-ggplot(aes(x= focus, y = meanB, fill = alcohol), data = satis_tib ) +
  facet_grid(~ satis, 
             labeller = labeller(satis = rating_labs))+
   geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=ci_lo, ymax=ci_hi), width=.2,
                 position=position_dodge(.9))+
  scale_y_continuous(name="satisfaction", 
                     limits=c(0, 100), 
                     breaks =   seq(0,100, 25)) +
  scale_x_discrete(name = "focus" )+
  scale_fill_manual(values = gender_pal) +
  theme_minimal()

ggplotly(satis_plot)

Expected and perceived ratings of satisfaction errorbars indicate 95% CI (+/- 1.96 SE)

# model m A --> M
msm <- lm(satisE ~ focus*alcohol, data = satis_dat)

# model A--> B
msy <-lm(satis ~ satisE + focus*alcohol,  data = satis_dat)

broom::tidy(msm)%>%
  kable(digits = 3, caption = "effect of alcoholand focus on expected satisfaction")%>%
  kable_styling(full_width = TRUE)
effect of alcoholand focus on expected satisfaction
term estimate std.error statistic p.value
(Intercept) 55.500 3.088 17.972 0.000
focushedonic -1.893 4.444 -0.426 0.671
focusutilitarian -0.203 4.155 -0.049 0.961
alcoholnon-alcoholic -8.500 4.367 -1.946 0.053
focushedonic:alcoholnon-alcoholic 5.955 6.183 0.963 0.337
focusutilitarian:alcoholnon-alcoholic 7.803 6.028 1.294 0.197
med_satis_focusU <- mediate(msm, msy, 
                           sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "satisE",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "utilitarian")

summary(med_satis_focusU)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME              0.682       -0.338         2.20    0.22
## ADE              -0.476       -6.678         5.83    0.88
## Total Effect      0.206       -6.527         6.45    0.94
## Prop. Mediated    3.312       -3.778         4.27    0.88
## 
## Sample Size Used: 187 
## 
## 
## Simulations: 1000
plot(med_satis_focusU)

med_satis_focusH <- mediate(msm, msy, 
                           sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "satisE",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "hedonic")

summary(med_satis_focusH)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME              0.194       -1.087         1.61    0.78
## ADE              -4.819      -12.295         2.51    0.22
## Total Effect     -4.625      -12.059         2.96    0.26
## Prop. Mediated   -0.042       -1.181         1.09    0.90
## 
## Sample Size Used: 187 
## 
## 
## Simulations: 1000
plot(med_satis_focusH)

med_satis_alcohol <- mediate(msm, msy, 
                             sims = 1000, 
                              boot.ci.type = "perc", 
                              treat = "alcohol", 
                              mediator = "satisE",
                              robustSE = TRUE, 
                              boot = TRUE,
                              conf.level = 0.95, 
                              covariates = "focus", 
                              control.value = "alcoholic", 
                              treat.value = "non-alcoholic")

summary(med_satis_alcohol)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME            -0.7113      -2.1322         0.17   0.150   
## ADE             -7.0777     -12.2471        -1.40   0.012 * 
## Total Effect    -7.7890     -13.0529        -2.07   0.008 **
## Prop. Mediated   0.0913      -0.0279         0.40   0.154   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 187 
## 
## 
## Simulations: 1000
plot(med_satis_alcohol)

Willingness to pay

  • no effect of focus on willingness to pay
  • effect of alcoholcontent on wtp ( pp would pay slightly more for alcoholic beer)
# model m A --> M
mwm <- lm(satis ~ focus*alcohol, data = satis_dat)

# model A--> B
mwy <-lm(wtp ~ satis + focus*alcohol,  data = satis_dat)


med_wtp_focusH <- mediate(mwm, mwy, 
                          sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "satis",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "hedonic")

summary(med_wtp_focusH)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME            -0.0943      -0.2411         0.05    0.21
## ADE             -0.0658      -0.3619         0.20    0.63
## Total Effect    -0.1601      -0.4736         0.15    0.31
## Prop. Mediated   0.5889      -3.6505         4.28    0.36
## 
## Sample Size Used: 187 
## 
## 
## Simulations: 1000
plot(med_wtp_focusH)

med_wtp_focusU <- mediate(mwm, mwy, 
                          sims = 1000, 
                            boot = TRUE,
                            boot.ci.type = "perc", 
                            treat = "focus", 
                            mediator = "satis",
                            covariates = "alcohol", 
                            robustSE = TRUE,  
                            conf.level = 0.95,
                            control.value = "control", 
                            treat.value = "utilitarian")

summary(med_wtp_focusU)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME             0.0042      -0.1337         0.14    0.95
## ADE             -0.1707      -0.4372         0.10    0.20
## Total Effect    -0.1665      -0.4779         0.15    0.28
## Prop. Mediated  -0.0252      -3.1410         3.37    0.85
## 
## Sample Size Used: 187 
## 
## 
## Simulations: 1000
plot(med_wtp_focusU)

med_wtp_alcohol <- mediate(mwm, mwy, 
                           sims = 1000, 
                              boot.ci.type = "perc", 
                              treat = "alcohol", 
                              mediator = "satis",
                              robustSE = TRUE, 
                              boot = TRUE,
                              conf.level = 0.95, 
                             covariates = "focus",
                              control.value = "alcoholic", 
                              treat.value = "non-alcoholic")

summary(med_wtp_alcohol)
## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
## (Inference Conditional on the Covariate Values Specified in `covariates')
## 
##                Estimate 95% CI Lower 95% CI Upper p-value   
## ACME            -0.1588      -0.2755        -0.04   0.004 **
## ADE             -0.0949      -0.3235         0.13   0.402   
## Total Effect    -0.2536      -0.5000        -0.01   0.036 * 
## Prop. Mediated   0.6260       0.1091         2.85   0.036 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 187 
## 
## 
## Simulations: 1000
plot(med_wtp_alcohol)

wtp_model <- robust::lmRob(wtp ~ focus * alcohol, data = satis_dat)
summary(wtp_model)
## 
## Call:
## robust::lmRob(formula = wtp ~ focus * alcohol, data = satis_dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1655 -0.7333  0.1019  0.8345  2.1912 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             3.1655     0.1790  17.683   <2e-16 ***
## focushedonic                           -0.2674     0.2554  -1.047   0.2966    
## focusutilitarian                       -0.3567     0.2403  -1.484   0.1394    
## alcoholnon-alcoholic                   -0.4321     0.2510  -1.721   0.0869 .  
## focushedonic:alcoholnon-alcoholic       0.1903     0.3539   0.538   0.5915    
## focusutilitarian:alcoholnon-alcoholic   0.3233     0.3474   0.931   0.3532    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9126 on 181 degrees of freedom
## Multiple R-Squared: 0.03504 
## 
## Test for Bias:
##             statistic p-value
## M-estimate     14.972 0.02048
## LS-estimate    -6.361 1.00000
wtp_tib <- satis_dat %>%
  group_by(focus, alcohol) %>%
  summarise(n = n(),
            mean_wtp = mean(wtp, na.rm = TRUE), 
            sd_wtp = sd(wtp, na.rm = TRUE), 
            ci_lo = mean_wtp - 1.96 * (sd_wtp/sqrt(n)),
            ci_high = mean_wtp + 1.96 *(sd_wtp/sqrt(n)))

wtp_tib %>%
  kable(digits = 2, caption = "Participants' mean willingness to pay in different alcohol and focus conditions.") %>%
  kable_styling(full_width = TRUE)
Participants’ mean willingness to pay in different alcohol and focus conditions.
focus alcohol n mean_wtp sd_wtp ci_lo ci_high
control alcoholic 30 3.13 0.90 2.81 3.46
control non-alcoholic 30 2.73 0.83 2.44 3.03
hedonic alcoholic 28 2.89 0.88 2.57 3.22
hedonic non-alcoholic 32 2.66 0.97 2.32 2.99
utilitarian alcoholic 37 2.84 0.73 2.60 3.07
utilitarian non-alcoholic 30 2.70 0.88 2.39 3.01
wtp_bar <-ggplot(aes(x=focus, 
           y = mean_wtp, 
           fill = alcohol), 
           data= wtp_tib)+
  geom_bar(stat="identity", position=position_dodge())+
  geom_errorbar(aes(ymin=ci_lo, ymax=ci_high), width=.2,
                 position=position_dodge(.9))+
  scale_y_continuous(name="willingness to pay (£)", 
                     limits=c(0, 5), 
                     breaks =   seq(0,5, 1)) +
  scale_x_discrete(name = "focus" )+
  scale_fill_manual(values = gender_pal) +
  theme_minimal()

wtp_bar
Participants' willingness to pay for alcoholic and non-alcoholic beer. Error bars indicate confidence intervals (+/- 1.96 SE)

Participants’ willingness to pay for alcoholic and non-alcoholic beer. Error bars indicate confidence intervals (+/- 1.96 SE)

Cognitive Performance

Participants’ cognitive performance was assessed using the inspection time task (10.1016/j.neuroimage.2004.03.047) which measures participants’information processing. Participants completed the task just before, 30 minutes and 60 minutes after consumption.

-pp performance improves as duration of stimulus increases - effect of alcohol, stimulus duration and stage - interactions: alcoholstimulus stimulusstage alcohol*stimulus*stage - no effect of focus

# these are two methods of analysis, should do the same thing, buttakes too long to run

# m_itt <- lmerTest::lmer(correct ~ focus * alcohol *stage *stimulus + (stage*stimulus| subject), data=itt_dat, REML=F)
# 
# summary(m_itt)

# contrasts

itt_afx <- afex::aov_4(correct ~ focus*alcohol*stage*stimulus  + (stimulus*stage|subject), data = itt_dat)
summary(itt_afx)
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                              Sum Sq num Df Error SS den Df    F value    Pr(>F)
## (Intercept)                  5831.3      1   38.734    178 26797.4015 < 2.2e-16
## focus                           0.6      2   38.734    178     1.4063  0.247768
## alcohol                         1.8      1   38.734    178     8.1225  0.004888
## focus:alcohol                   0.7      2   38.734    178     1.5975  0.205284
## stimulus                       34.3     12   15.356   2136   397.8339 < 2.2e-16
## focus:stimulus                  0.3     24   15.356   2136     1.4808  0.062289
## alcohol:stimulus                0.2     12   15.356   2136     1.8983  0.030276
## focus:alcohol:stimulus          0.1     24   15.356   2136     0.5317  0.969485
## stage                           0.5      2    8.437    356     9.7672 7.421e-05
## focus:stage                     0.1      4    8.437    356     1.0657  0.373312
## alcohol:stage                   0.0      2    8.437    356     0.2461  0.781995
## focus:alcohol:stage             0.0      4    8.437    356     0.0856  0.986846
## stimulus:stage                  0.5     24   18.452   4272     4.8156 9.804e-14
## focus:stimulus:stage            0.2     48   18.452   4272     0.8259  0.798308
## alcohol:stimulus:stage          0.2     24   18.452   4272     1.5525  0.041858
## focus:alcohol:stimulus:stage    0.2     48   18.452   4272     0.7592  0.887985
##                                 
## (Intercept)                  ***
## focus                           
## alcohol                      ** 
## focus:alcohol                   
## stimulus                     ***
## focus:stimulus               .  
## alcohol:stimulus             *  
## focus:alcohol:stimulus          
## stage                        ***
## focus:stage                     
## alcohol:stage                   
## focus:alcohol:stage             
## stimulus:stage               ***
## focus:stimulus:stage            
## alcohol:stimulus:stage       *  
## focus:alcohol:stimulus:stage    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                              Test statistic    p-value
## stimulus                            0.00230 0.0000e+00
## focus:stimulus                      0.00230 0.0000e+00
## alcohol:stimulus                    0.00230 0.0000e+00
## focus:alcohol:stimulus              0.00230 0.0000e+00
## stage                               0.58653 3.1174e-21
## focus:stage                         0.58653 3.1174e-21
## alcohol:stage                       0.58653 3.1174e-21
## focus:alcohol:stage                 0.58653 3.1174e-21
## stimulus:stage                      0.00023 0.0000e+00
## focus:stimulus:stage                0.00023 0.0000e+00
## alcohol:stimulus:stage              0.00023 0.0000e+00
## focus:alcohol:stimulus:stage        0.00023 0.0000e+00
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                               GG eps Pr(>F[GG])    
## stimulus                     0.43240  < 2.2e-16 ***
## focus:stimulus               0.43240  0.1380588    
## alcohol:stimulus             0.43240  0.0893189 .  
## focus:alcohol:stimulus       0.43240  0.8740470    
## stage                        0.70748  0.0005176 ***
## focus:stage                  0.70748  0.3619277    
## alcohol:stage                0.70748  0.7031497    
## focus:alcohol:stage          0.70748  0.9624416    
## stimulus:stage               0.52543  3.705e-08 ***
## focus:stimulus:stage         0.52543  0.7126193    
## alcohol:stimulus:stage       0.52543  0.0944288 .  
## focus:alcohol:stimulus:stage 0.52543  0.7985395    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                                 HF eps    Pr(>F[HF])
## stimulus                     0.4469239 4.817097e-240
## focus:stimulus               0.4469239  1.351004e-01
## alcohol:stimulus             0.4469239  8.677605e-02
## focus:alcohol:stimulus       0.4469239  8.789647e-01
## stage                        0.7114920  5.039157e-04
## focus:stage                  0.7114920  3.621441e-01
## alcohol:stage                0.7114920  7.044778e-01
## focus:alcohol:stage          0.7114920  9.629908e-01
## stimulus:stage               0.5681699  1.156033e-08
## focus:stimulus:stage         0.5681699  7.226680e-01
## alcohol:stimulus:stage       0.5681699  8.748527e-02
## focus:alcohol:stimulus:stage 0.5681699  8.098047e-01
# need to take effects apart and set contrasts

# the main effect of alcohol

alcohol_emm <- emmeans::emmeans(itt_afx, ~alcohol, model = "multivariate")
alcohol_emm # shows us the means
##  alcohol       emmean      SE  df lower.CL upper.CL
##  alcoholic      0.887 0.00777 178    0.872    0.903
##  non-alcoholic  0.919 0.00783 178    0.903    0.934
## 
## Results are averaged over the levels of: focus, stage, stimulus 
## Confidence level used: 0.95
# the main effect of stage
stage_emm <- emmeans::emmeans(itt_afx, ~stage, model = "multivariate")
stage_emm
##  stage   emmean      SE  df lower.CL upper.CL
##  stage.1  0.909 0.00632 178    0.896    0.921
##  stage.2  0.908 0.00553 178    0.897    0.919
##  stage.3  0.892 0.00637 178    0.879    0.904
## 
## Results are averaged over the levels of: focus, alcohol, stimulus 
## Confidence level used: 0.95
emmeans::contrast(stage_emm, method = "trt.vs.ctrl", ref = 1, adjust = "holm")
##  contrast           estimate      SE  df t.ratio p.value
##  stage.2 - stage.1 -0.000635 0.00473 178 -0.134  0.8933 
##  stage.3 - stage.1 -0.017376 0.00543 178 -3.198  0.0033 
## 
## Results are averaged over the levels of: focus, alcohol, stimulus 
## P value adjustment: holm method for 2 tests
#main effect of stimulus
stimulus_emm <- emmeans::emmeans(itt_afx, ~stimulus, model = "multivariate")
stimulus_emm
##  stimulus emmean      SE  df lower.CL upper.CL
##  stim1     0.701 0.00747 178    0.686    0.716
##  stim2     0.846 0.00784 178    0.830    0.861
##  stim3     0.857 0.00733 178    0.843    0.871
##  stim4     0.862 0.00770 178    0.846    0.877
##  stim5     0.921 0.00682 178    0.908    0.934
##  stim6     0.924 0.00623 178    0.911    0.936
##  stim7     0.934 0.00610 178    0.922    0.946
##  stim8     0.944 0.00574 178    0.933    0.955
##  stim9     0.945 0.00586 178    0.934    0.957
##  stim10    0.948 0.00586 178    0.936    0.959
##  stim11    0.947 0.00597 178    0.935    0.959
##  stim12    0.954 0.00550 178    0.943    0.965
##  stim13    0.956 0.00566 178    0.945    0.967
## 
## Results are averaged over the levels of: focus, alcohol, stage 
## Confidence level used: 0.95
emmeans::contrast(stimulus_emm, method = "poly", adjust = "holm")
##  contrast  estimate     SE  df t.ratio p.value
##  linear        2.76 0.0932 178  29.566 <.0001 
##  quadratic    -5.21 0.1856 178 -28.062 <.0001 
##  cubic         1.32 0.0952 178  13.873 <.0001 
##  quartic      -7.12 0.9739 178  -7.307 <.0001 
##  degree 5      2.40 0.2484 178   9.660 <.0001 
##  degree 6     -3.96 0.3900 178 -10.149 <.0001 
## 
## Results are averaged over the levels of: focus, alcohol, stage 
## P value adjustment: holm method for 6 tests
# interaction stimulus:stage

inter_emm <- emmeans::emmeans(itt_afx, c("stage", "stimulus"), model = "multivariate")
inter_emm
##  stage   stimulus emmean      SE  df lower.CL upper.CL
##  stage.1 stim1     0.672 0.01041 178    0.652    0.693
##  stage.2 stim1     0.716 0.01077 178    0.695    0.738
##  stage.3 stim1     0.715 0.01037 178    0.694    0.735
##  stage.1 stim2     0.863 0.00923 178    0.844    0.881
##  stage.2 stim2     0.847 0.00921 178    0.829    0.865
##  stage.3 stim2     0.828 0.00976 178    0.809    0.847
##  stage.1 stim3     0.865 0.00876 178    0.848    0.883
##  stage.2 stim3     0.870 0.00889 178    0.852    0.887
##  stage.3 stim3     0.836 0.00896 178    0.818    0.854
##  stage.1 stim4     0.881 0.00909 178    0.863    0.899
##  stage.2 stim4     0.855 0.00907 178    0.837    0.873
##  stage.3 stim4     0.849 0.00973 178    0.830    0.868
##  stage.1 stim5     0.924 0.00835 178    0.907    0.940
##  stage.2 stim5     0.928 0.00749 178    0.914    0.943
##  stage.3 stim5     0.911 0.00841 178    0.894    0.927
##  stage.1 stim6     0.936 0.00682 178    0.922    0.949
##  stage.2 stim6     0.928 0.00727 178    0.914    0.943
##  stage.3 stim6     0.907 0.00783 178    0.891    0.922
##  stage.1 stim7     0.941 0.00780 178    0.925    0.956
##  stage.2 stim7     0.937 0.00638 178    0.924    0.949
##  stage.3 stim7     0.926 0.00749 178    0.911    0.940
##  stage.1 stim8     0.951 0.00750 178    0.936    0.966
##  stage.2 stim8     0.948 0.00575 178    0.937    0.959
##  stage.3 stim8     0.932 0.00814 178    0.916    0.948
##  stage.1 stim9     0.951 0.00758 178    0.936    0.966
##  stage.2 stim9     0.953 0.00625 178    0.941    0.965
##  stage.3 stim9     0.932 0.00689 178    0.919    0.946
##  stage.1 stim10    0.956 0.00723 178    0.941    0.970
##  stage.2 stim10    0.955 0.00653 178    0.942    0.968
##  stage.3 stim10    0.932 0.00727 178    0.918    0.947
##  stage.1 stim11    0.952 0.00783 178    0.937    0.968
##  stage.2 stim11    0.953 0.00610 178    0.941    0.965
##  stage.3 stim11    0.936 0.00765 178    0.921    0.952
##  stage.1 stim12    0.963 0.00717 178    0.949    0.977
##  stage.2 stim12    0.955 0.00566 178    0.944    0.966
##  stage.3 stim12    0.943 0.00713 178    0.929    0.957
##  stage.1 stim13    0.962 0.00665 178    0.949    0.975
##  stage.2 stim13    0.962 0.00626 178    0.950    0.974
##  stage.3 stim13    0.944 0.00715 178    0.930    0.958
## 
## Results are averaged over the levels of: focus, alcohol 
## Confidence level used: 0.95
emmeans::contrast(
  inter_emm,
  c(stage = "trt.vs.ctrl", stimulus = "poly"),
  rel = 1,
  adjust = "holm"
  )
##  contrast                       estimate     SE  df t.ratio p.value
##  stage.2 stim1 - stage.1 stim1    0.0439 0.0135 178  3.257  0.0027 
##  stage.3 stim1 - stage.1 stim1    0.0421 0.0137 178  3.079  0.0027 
##  stage.1 stim2 - stage.1 stim1    0.1902 0.0106 178 18.003  <.0001 
##  stage.2 stim2 - stage.1 stim1    0.1742 0.0116 178 15.007  <.0001 
##  stage.3 stim2 - stage.1 stim1    0.1557 0.0117 178 13.329  <.0001 
##  stage.1 stim3 - stage.1 stim1    0.1930 0.0101 178 19.207  <.0001 
##  stage.2 stim3 - stage.1 stim1    0.1973 0.0104 178 18.880  <.0001 
##  stage.3 stim3 - stage.1 stim1    0.1634 0.0111 178 14.660  <.0001 
##  stage.1 stim4 - stage.1 stim1    0.2082 0.0110 178 18.840  <.0001 
##  stage.2 stim4 - stage.1 stim1    0.1829 0.0113 178 16.166  <.0001 
##  stage.3 stim4 - stage.1 stim1    0.1763 0.0119 178 14.801  <.0001 
##  stage.1 stim5 - stage.1 stim1    0.2515 0.0110 178 22.820  <.0001 
##  stage.2 stim5 - stage.1 stim1    0.2560 0.0112 178 22.949  <.0001 
##  stage.3 stim5 - stage.1 stim1    0.2382 0.0109 178 21.847  <.0001 
##  stage.1 stim6 - stage.1 stim1    0.2634 0.0106 178 24.888  <.0001 
##  stage.2 stim6 - stage.1 stim1    0.2560 0.0113 178 22.664  <.0001 
##  stage.3 stim6 - stage.1 stim1    0.2344 0.0114 178 20.518  <.0001 
##  stage.1 stim7 - stage.1 stim1    0.2682 0.0107 178 24.987  <.0001 
##  stage.2 stim7 - stage.1 stim1    0.2641 0.0105 178 25.037  <.0001 
##  stage.3 stim7 - stage.1 stim1    0.2531 0.0112 178 22.671  <.0001 
##  stage.1 stim8 - stage.1 stim1    0.2786 0.0108 178 25.780  <.0001 
##  stage.2 stim8 - stage.1 stim1    0.2756 0.0102 178 26.988  <.0001 
##  stage.3 stim8 - stage.1 stim1    0.2599 0.0116 178 22.321  <.0001 
##  stage.1 stim9 - stage.1 stim1    0.2784 0.0110 178 25.341  <.0001 
##  stage.2 stim9 - stage.1 stim1    0.2806 0.0111 178 25.329  <.0001 
##  stage.3 stim9 - stage.1 stim1    0.2598 0.0111 178 23.461  <.0001 
##  stage.1 stim10 - stage.1 stim1   0.2832 0.0110 178 25.756  <.0001 
##  stage.2 stim10 - stage.1 stim1   0.2827 0.0106 178 26.659  <.0001 
##  stage.3 stim10 - stage.1 stim1   0.2598 0.0117 178 22.153  <.0001 
##  stage.1 stim11 - stage.1 stim1   0.2798 0.0114 178 24.602  <.0001 
##  stage.2 stim11 - stage.1 stim1   0.2806 0.0108 178 25.919  <.0001 
##  stage.3 stim11 - stage.1 stim1   0.2640 0.0117 178 22.576  <.0001 
##  stage.1 stim12 - stage.1 stim1   0.2907 0.0109 178 26.560  <.0001 
##  stage.2 stim12 - stage.1 stim1   0.2828 0.0112 178 25.360  <.0001 
##  stage.3 stim12 - stage.1 stim1   0.2706 0.0117 178 23.037  <.0001 
##  stage.1 stim13 - stage.1 stim1   0.2892 0.0108 178 26.842  <.0001 
##  stage.2 stim13 - stage.1 stim1   0.2896 0.0109 178 26.655  <.0001 
##  stage.3 stim13 - stage.1 stim1   0.2714 0.0115 178 23.677  <.0001 
## 
## Results are averaged over the levels of: focus, alcohol 
## P value adjustment: holm method for 38 tests
itt_dat$stimulus <- as.factor(itt_dat$stimulus) %>%
  fct_relevel(.,c("stim1", "stim2", "stim3", "stim4", "stim5", "stim6", "stim7", "stim8", "stim9", "stim10", "stim11", "stim12","stim13"))


itt_plot<-ggplot(data = itt_dat, aes(x = stimulus, y = correct, color = alcohol))+
  facet_grid(rows = vars(focus),
             cols = vars(stage))+
  stat_summary(fun=mean, geom="point") +
  scale_y_continuous(name="% correct", limits=c(0.5, 1.0), breaks = seq(0.5,1,0.15)) +
  scale_x_discrete(name = "stimulus duration\n(ms)", labels =c("19", "25", "31", "37", "44", "50", "62", "75", "87", "100", "125", "150", "200") )+ 
  scale_color_manual(values = gender_pal) +
  theme(axis.text.x = element_text( color="black", 
                           size=2, angle=45))+
  theme_minimal()

ggplotly(itt_plot, height = 500, width=1200)

Participants’ cognitive performance

Mood

mood_data<- mood_dat %>%
  dplyr::select(subject, stage, calm, content, alert)%>%
  pivot_longer(cols = c(calm, content, alert),
               names_to = "factor",
               values_to = "rating")%>%
  mutate(condition = subject%/% 1000)%>%
  mutate(alcohol = if_else((condition%%2)==1, "alcoholic", "non-alcoholic"))%>%
  mutate(focus = case_when(condition == 1 ~ 'control',
                           condition == 2 ~ 'control',
                           condition == 3 ~ 'hedonic',
                           condition == 4 ~ 'hedonic',
                           condition == 5 ~ 'utilitarian',
                           condition == 6 ~ 'utilitarian'))

ggplot(data = mood_data, aes(x = stage, y = rating, color = alcohol))+
    facet_grid(rows = vars(focus),
             cols = vars(factor))+
  stat_summary(fun=mean, geom="point") +
  scale_y_continuous(name="rating", limits=c(35, 48), breaks = seq(0,50,5)) +
  scale_x_discrete(name = "stage", labels = c("1", "2", "3") )+ 
  scale_color_manual(values = gender_pal) +
  theme(axis.text.x = element_text( color="black", 
                           size=2, angle=45))+
  theme_light()
Participants'mood accross the stages and conditions

Participants’mood accross the stages and conditions

mood_data$focus <- as.factor(mood_data$focus)

# afex mixed effects model
mood_afx <- afex::aov_4(rating ~ focus*alcohol*stage*factor  + (factor*stage|subject), data = mood_data)
summary(mood_afx)
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                             Sum Sq num Df Error SS den Df   F value    Pr(>F)
## (Intercept)                2505627      1   183226    179 2447.8317 < 2.2e-16
## focus                         1322      2   183226    179    0.6456 0.5255832
## alcohol                       1250      1   183226    179    1.2217 0.2705194
## focus:alcohol                 9356      2   183226    179    4.5701 0.0115939
## factor                       33302      2    72383    358   82.3551 < 2.2e-16
## focus:factor                  1339      4    72383    358    1.6552 0.1599128
## alcohol:factor                 888      2    72383    358    2.1963 0.1127118
## focus:alcohol:factor           683      4    72383    358    0.8446 0.4975886
## stage                         1456      2    39514    358    6.5970 0.0015363
## focus:stage                    576      4    39514    358    1.3050 0.2676981
## alcohol:stage                  416      2    39514    358    1.8827 0.1536885
## focus:alcohol:stage            402      4    39514    358    0.9107 0.4577106
## factor:stage                  8530      4    42998    716   35.5088 < 2.2e-16
## focus:factor:stage             182      8    42998    716    0.3780 0.9324598
## alcohol:factor:stage          1201      4    42998    716    4.9981 0.0005606
## focus:alcohol:factor:stage     477      8    42998    716    0.9932 0.4398359
##                               
## (Intercept)                ***
## focus                         
## alcohol                       
## focus:alcohol              *  
## factor                     ***
## focus:factor                  
## alcohol:factor                
## focus:alcohol:factor          
## stage                      ** 
## focus:stage                   
## alcohol:stage                 
## focus:alcohol:stage           
## factor:stage               ***
## focus:factor:stage            
## alcohol:factor:stage       ***
## focus:alcohol:factor:stage    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##                            Test statistic    p-value
## factor                            0.78985 0.00000000
## focus:factor                      0.78985 0.00000000
## alcohol:factor                    0.78985 0.00000000
## focus:alcohol:factor              0.78985 0.00000000
## stage                             0.92481 0.00095205
## focus:stage                       0.92481 0.00095205
## alcohol:stage                     0.92481 0.00095205
## focus:alcohol:stage               0.92481 0.00095205
## factor:stage                      0.45122 0.00000000
## focus:factor:stage                0.45122 0.00000000
## alcohol:factor:stage              0.45122 0.00000000
## focus:alcohol:factor:stage        0.45122 0.00000000
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                             GG eps Pr(>F[GG])    
## factor                     0.82634  < 2.2e-16 ***
## focus:factor               0.82634   0.171718    
## alcohol:factor             0.82634   0.122566    
## focus:alcohol:factor       0.82634   0.479646    
## stage                      0.93007   0.002021 ** 
## focus:stage                0.93007   0.269465    
## alcohol:stage              0.93007   0.156820    
## focus:alcohol:stage        0.93007   0.452394    
## factor:stage               0.75490  < 2.2e-16 ***
## focus:factor:stage         0.75490   0.894105    
## alcohol:factor:stage       0.75490   0.001944 ** 
## focus:alcohol:factor:stage 0.75490   0.429238    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                               HF eps   Pr(>F[HF])
## factor                     0.8330648 1.429346e-25
## focus:factor               0.8330648 1.712519e-01
## alcohol:factor             0.8330648 1.221801e-01
## focus:alcohol:factor       0.8330648 4.804104e-01
## stage                      0.9394397 1.947862e-03
## focus:stage                0.9394397 2.692367e-01
## alcohol:stage              0.9394397 1.564039e-01
## focus:alcohol:stage        0.9394397 4.531321e-01
## factor:stage               0.7693046 2.017256e-21
## focus:factor:stage         0.7693046 8.969476e-01
## alcohol:factor:stage       0.7693046 1.806356e-03
## focus:alcohol:factor:stage 0.7693046 4.299908e-01
# main effect stage

stage_emm <- emmeans::emmeans(mood_afx, "stage", model = "multivariate")
stage_emm
##  stage     emmean    SE  df lower.CL upper.CL
##  response1   37.7 0.905 179     35.9     39.4
##  response2   39.0 0.802 179     37.5     40.6
##  response3   39.9 0.888 179     38.2     41.7
## 
## Results are averaged over the levels of: focus, alcohol, factor 
## Confidence level used: 0.95
emmeans::contrast(
  stage_emm,
  stage = "trt.vs.ctrl",
  ref = 1,
  adjust = "holm"
  )
##  contrast         estimate    SE  df t.ratio p.value
##  response1 effect   -1.222 0.412 179 -2.967  0.0074 
##  response2 effect    0.165 0.334 179  0.494  0.6216 
##  response3 effect    1.056 0.344 179  3.073  0.0074 
## 
## Results are averaged over the levels of: focus, alcohol, factor 
## P value adjustment: holm method for 3 tests
# main effect mood type = factor
fact_emm <- emmeans::emmeans(mood_afx, "factor", model = "multivariate")
fact_emm
##  factor  emmean    SE  df lower.CL upper.CL
##  calm      34.5 0.868 179     32.8     36.2
##  content   37.1 0.995 179     35.1     39.0
##  alert     45.0 0.917 179     43.2     46.9
## 
## Results are averaged over the levels of: focus, alcohol, stage 
## Confidence level used: 0.95
emmeans::contrast(
  fact_emm,
  stage = "pairwise",
  adjust = "holm"
  )
##  contrast       estimate    SE  df t.ratio p.value
##  calm effect       -4.35 0.581 179 -7.484  <.0001 
##  content effect    -1.82 0.379 179 -4.804  <.0001 
##  alert effect       6.17 0.501 179 12.324  <.0001 
## 
## Results are averaged over the levels of: focus, alcohol, stage 
## P value adjustment: holm method for 3 tests
# interaction focus:alcohol

foc_alc_emm <- emmeans::emmeans(mood_afx, c("focus", "alcohol"), model = "multivariate")
foc_alc_emm
##  focus       alcohol       emmean   SE  df lower.CL upper.CL
##  control     alcoholic       37.2 1.92 179     33.5     41.0
##  hedonic     alcoholic       43.5 2.02 179     39.6     47.5
##  utilitarian alcoholic       38.5 1.80 179     34.9     42.0
##  control     non-alcoholic   38.0 1.95 179     34.2     41.9
##  hedonic     non-alcoholic   35.1 1.92 179     31.3     38.9
##  utilitarian non-alcoholic   40.9 1.95 179     37.1     44.8
## 
## Results are averaged over the levels of: stage, factor 
## Confidence level used: 0.95
emmeans::contrast(
  foc_alc_emm,
  stage = "poly",
  adjust = "holm"
  )
##  contrast                           estimate   SE  df t.ratio p.value
##  control alcoholic effect             -1.637 1.75 179 -0.935  1.0000 
##  hedonic alcoholic effect              4.648 1.82 179  2.549  0.0699 
##  utilitarian alcoholic effect         -0.405 1.67 179 -0.243  1.0000 
##  (control non-alcoholic) effect       -0.840 1.77 179 -0.474  1.0000 
##  (hedonic non-alcoholic) effect       -3.795 1.75 179 -2.168  0.1573 
##  (utilitarian non-alcoholic) effect    2.029 1.77 179  1.144  1.0000 
## 
## Results are averaged over the levels of: stage, factor 
## P value adjustment: holm method for 6 tests
#interaction factor:stage
fact_stage_emm <- emmeans::emmeans(mood_afx, c("factor", "stage"), model = "multivariate")
fact_stage_emm
##  factor  stage     emmean    SE  df lower.CL upper.CL
##  calm    response1   36.7 1.110 179     34.5     38.9
##  content response1   36.7 1.150 179     34.4     38.9
##  alert   response1   39.6 1.054 179     37.5     41.7
##  calm    response2   33.3 0.996 179     31.3     35.2
##  content response2   36.3 1.030 179     34.3     38.4
##  alert   response2   47.5 0.988 179     45.6     49.5
##  calm    response3   33.6 1.128 179     31.4     35.8
##  content response3   38.2 1.076 179     36.1     40.3
##  alert   response3   48.0 1.064 179     45.9     50.1
## 
## Results are averaged over the levels of: focus, alcohol 
## Confidence level used: 0.95
emmeans::contrast(
  fact_stage_emm,
  interaction = c(stage = "pairwise", fact = "trt.vs.ctrl"),
  ref = 1,
  adjust = "holm"
  )
##  stage_pairwise        factor_trt.vs.ctrl estimate    SE  df t.ratio p.value
##  response1 - response2 content - calm       -3.151 1.219 179 -2.584  0.0317 
##  response1 - response3 content - calm       -4.664 1.246 179 -3.742  0.0010 
##  response2 - response3 content - calm       -1.513 0.996 179 -1.520  0.2608 
##  response1 - response2 alert - calm        -11.415 1.374 179 -8.311  <.0001 
##  response1 - response3 alert - calm        -11.593 1.472 179 -7.875  <.0001 
##  response2 - response3 alert - calm         -0.178 1.232 179 -0.144  0.8853 
## 
## Results are averaged over the levels of: focus, alcohol 
## P value adjustment: holm method for 6 tests
#interaction alcohol:factor:stage

fact_stage_alc_emm <- emmeans::emmeans(mood_afx, c("factor", "alcohol", "stage"), model = "multivariate")
fact_stage_alc_emm
##  factor  alcohol       stage     emmean   SE  df lower.CL upper.CL
##  calm    alcoholic     response1   38.7 1.56 179     35.6     41.8
##  content alcoholic     response1   36.4 1.62 179     33.2     39.6
##  alert   alcoholic     response1   38.7 1.48 179     35.8     41.6
##  calm    non-alcoholic response1   34.8 1.58 179     31.6     37.9
##  content non-alcoholic response1   36.9 1.64 179     33.7     40.1
##  alert   non-alcoholic response1   40.5 1.50 179     37.5     43.4
##  calm    alcoholic     response2   34.3 1.40 179     31.6     37.1
##  content alcoholic     response2   35.3 1.45 179     32.5     38.2
##  alert   alcoholic     response2   49.9 1.39 179     47.2     52.7
##  calm    non-alcoholic response2   32.2 1.42 179     29.4     35.0
##  content non-alcoholic response2   37.3 1.47 179     34.4     40.2
##  alert   non-alcoholic response2   45.1 1.41 179     42.4     47.9
##  calm    alcoholic     response3   35.1 1.59 179     32.0     38.3
##  content alcoholic     response3   38.9 1.51 179     36.0     41.9
##  alert   alcoholic     response3   50.2 1.50 179     47.3     53.2
##  calm    non-alcoholic response3   32.0 1.61 179     28.9     35.2
##  content non-alcoholic response3   37.4 1.53 179     34.4     40.4
##  alert   non-alcoholic response3   45.8 1.51 179     42.9     48.8
## 
## Results are averaged over the levels of: focus 
## Confidence level used: 0.95
emmeans::contrast(
  fact_stage_alc_emm,
  interaction = c(stage = "pairwise", fact = "trt.vs.ctrl"),
  ref = 1,
  adjust = "holm"
  )
##  stage_pairwise        factor_trt.vs.ctrl alcohol_pairwise            estimate
##  response1 - response2 content - calm     alcoholic - (non-alcoholic)   -0.346
##  response1 - response3 content - calm     alcoholic - (non-alcoholic)   -2.872
##  response2 - response3 content - calm     alcoholic - (non-alcoholic)   -2.526
##  response1 - response2 alert - calm       alcoholic - (non-alcoholic)   -8.413
##  response1 - response3 alert - calm       alcoholic - (non-alcoholic)   -7.044
##  response2 - response3 alert - calm       alcoholic - (non-alcoholic)    1.369
##    SE  df t.ratio p.value
##  2.44 179 -0.142  1.0000 
##  2.49 179 -1.152  0.8256 
##  1.99 179 -1.268  0.8256 
##  2.75 179 -3.063  0.0152 
##  2.94 179 -2.392  0.0889 
##  2.46 179  0.556  1.0000 
## 
## Results are averaged over the levels of: focus 
## P value adjustment: holm method for 6 tests